# Load necessary librarieslibrary(forecast)library(tseries)library(tidyverse)library(tidymodels)library(kableExtra)library(modeltime)library(timetk)library(cowplot)# setting the seedset.seed(123)# Loading Data ----------------------------------------------------------------daily_deaths <-read_csv("data/new_daily_deaths.csv") %>%select(daily_deaths)east <-read_csv("data/east_daily.csv")midwest <-read_csv("data/midwest_daily.csv")south <-read_csv("data/south_daily.csv")west <-read_csv("data/west_daily.csv")
2. Data Splitting
Code
# Data Splitting ---------------------------------------------------------------east_splits <-time_series_split(east, initial ="1022 days", assess ="114 days")east_train <-training(east_splits)east_test <-testing(east_splits)midwest_splits <-time_series_split(midwest, initial ="1022 days", assess ="114 days")midwest_train <-training(midwest_splits)midwest_test <-testing(midwest_splits)south_splits <-time_series_split(south, initial ="1022 days", assess ="114 days")south_train <-training(south_splits)south_test <-testing(south_splits)west_splits <-time_series_split(west, initial ="1022 days", assess ="114 days")west_train <-training(west_splits)west_test <-testing(west_splits)## Convert Into Time Series ----------------------------------------------------east_train_ts <-ts(east_train %>%select(daily_deaths))east_test_ts <-ts(east_test %>%select(daily_deaths))east_ts <-ts(east %>%select(daily_deaths))midwest_train_ts <-ts(midwest_train %>%select(daily_deaths))midwest_test_ts <-ts(midwest_test %>%select(daily_deaths))midwest_ts <-ts(midwest %>%select(daily_deaths))south_train_ts <-ts(south_train %>%select(daily_deaths))south_test_ts <-ts(south_test %>%select(daily_deaths))south_ts <-ts(south %>%select(daily_deaths))west_train_ts <-ts(west_train %>%select(daily_deaths))west_test_ts <-ts(west_test %>%select(daily_deaths))west_ts <-ts(west %>%select(daily_deaths))## Plotting Training and Testing -----------------------------------------------par(mfrow=c(2,2))ggplot() +geom_line(data = east_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = east_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = east_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = east_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="East Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = midwest_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = midwest_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = midwest_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = midwest_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="Midwest Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = south_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = south_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = south_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = south_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="South Region Training vs. Testing Data")
Code
ggplot() +geom_line(data = west_test, aes(x = date, y = daily_deaths, color ="turquoise")) +geom_point(data = west_test, aes(x = date, y = daily_deaths, color ="turquoise"), size = .5) +geom_line(data = west_train, aes(x = date, y = daily_deaths, color ="orange")) +geom_point(data = west_train, aes(x = date, y = daily_deaths, color ="orange"), size = .5) +scale_color_manual(name="Data Set", values =c("turquoise", "orange"), labels =c("Training", "Testing")) +labs(y ="Daily Deaths",x ="Date",title ="West Region Training vs. Testing Data")
3. Testing for Stationarity
Code
### Using Augmented Dickey-Fuller Test -----------------------------------------adf_east <-adf.test(east_ts)print(adf_east) # p-value is 0.01 < 0.05, implying it is stationary
Augmented Dickey-Fuller Test
data: east_ts
Dickey-Fuller = -4.5566, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_midwest <-adf.test(midwest_ts)print(adf_midwest) # p-value is 0.4919 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: midwest_ts
Dickey-Fuller = -2.5984, Lag order = 10, p-value = 0.325
alternative hypothesis: stationary
Code
midwest_ts <-diff(midwest_ts)adf_midwest <-adf.test(midwest_ts)print(adf_midwest) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: midwest_ts
Dickey-Fuller = -13.43, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_south <-adf.test(south_ts)print(adf_south) # p-value is 0.3018 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: south_ts
Dickey-Fuller = -2.6533, Lag order = 10, p-value = 0.3018
alternative hypothesis: stationary
Code
south_ts <-diff(south_ts)adf_south <-adf.test(south_ts)print(adf_south) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: south_ts
Dickey-Fuller = -11.36, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
Code
adf_west <-adf.test(west_ts)print(adf_west) # p-value is 0.4098 > 0.05, implying it is NOT stationary
Augmented Dickey-Fuller Test
data: west_ts
Dickey-Fuller = -2.3981, Lag order = 10, p-value = 0.4098
alternative hypothesis: stationary
Code
west_ts <-diff(west_ts)adf_west <-adf.test(west_ts)print(adf_west) # p-value is 0.01 < 0.05, implying it is stationary with 1 diff
Augmented Dickey-Fuller Test
data: west_ts
Dickey-Fuller = -12.247, Lag order = 10, p-value = 0.01
alternative hypothesis: stationary
## ACF & PACF Plots ------------------------------------------------------------par(mfrow=c(3, 1))acf(east_ts, main ="East")pacf(east_ts, main ="East")plot(east_ts, main ="Daily Deaths: East")
Code
acf(midwest_ts, main ="Midwest")pacf(midwest_ts, main ="Midwest")plot(midwest_ts, main ="Daily Deaths: Midwest")
Code
acf(south_ts, main ="South")pacf(south_ts, main ="South")plot(south_ts, main ="Daily Deaths: South")
Code
acf(west_ts, main ="West")pacf(west_ts, main ="West")plot(south_ts, main ="Daily Deaths: South")
Code
plot(west_ts, main ="Daily Deaths: West")
5. Time Series Decomposition
Code
# total daily death decompositiondaily_ts <-ts(daily_deaths, frequency =365)decomposed_ts <-decompose(daily_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(daily_ts, main ="Total Original Time Series", ylab ="Value", col ="blue")plot(decomposed_ts$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_ts$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_ts$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
# East Decomposition -----------------------------------------------------------east_daily <- east %>%select(daily_deaths)east_ts <-ts(east_daily, frequency =365)decomposed_east <-decompose(east_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(east_ts, main ="East Original Time Series", ylab ="Value", col ="blue")plot(decomposed_east$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_east$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_east$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## Midwest Decomposition ----------------------------------------------------------midwest_daily <- midwest %>%select(daily_deaths)midwest_ts <-ts(midwest_daily, frequency =365)decomposed_midwest <-decompose(midwest_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(midwest_ts, main ="Midwest Original Time Series", ylab ="Value", col ="blue")plot(decomposed_midwest$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_midwest$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_midwest$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## South Decomposition ----------------------------------------------------------south_daily <- south %>%select(daily_deaths)south_ts <-ts(south_daily, frequency =365)decomposed_south <-decompose(south_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(south_ts, main ="South Original Time Series", ylab ="Value", col ="blue")plot(decomposed_south$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_south$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_south$random, main ="Residual Component", ylab ="Value", col ="purple")
Code
## West Decomposition ----------------------------------------------------------west_daily <- west %>%select(daily_deaths)west_ts <-ts(west_daily, frequency =365)decomposed_west <-decompose(west_ts)# Plot the original time series and the decomposed componentspar(mfrow=c(2,1))plot(west_ts, main ="West Original Time Series", ylab ="Value", col ="blue")plot(decomposed_west$trend, main ="Trend Component", ylab ="Value", col ="red")
Code
plot(decomposed_west$seasonal, main ="Seasonal Component", ylab ="Value", col ="green")plot(decomposed_west$random, main ="Residual Component", ylab ="Value", col ="purple")
6. ARIMA Model
6.1 Initial Model
Code
# Building ARIMA Model ---------------------------------------------------------east_arima <-arima(east_train_ts, order=c(1,0,1))midwest_arima <-arima(midwest_train_ts, order=c(1,1,1))south_arima <-arima(south_train_ts, order=c(1,1,1))west_arima <-arima(west_train_ts, order=c(1,1,1))## Summary & Forecast ---------------------------------------------------------summary(east_arima)
Call:
arima(x = east_train_ts, order = c(1, 0, 1))
Coefficients:
ar1 ma1 intercept
0.9176 -0.0388 197.3738
s.e. 0.0147 0.0493 39.7998
sigma^2 estimated as 12153: log likelihood = -6257.17, aic = 12522.33
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1151898 110.2407 64.66739 -Inf Inf 0.9918341 0.008515868
Code
summary(midwest_arima)
Call:
arima(x = midwest_train_ts, order = c(1, 1, 1))
Coefficients:
ar1 ma1
-0.3477 -1.0000
s.e. 0.0293 0.0026
sigma^2 estimated as 50173: log likelihood = -6970.99, aic = 13947.98
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.8748943 223.8834 138.9352 NaN Inf 0.6215031 -0.06581815
Code
summary(south_arima)
Call:
arima(x = south_train_ts, order = c(1, 1, 1))
Coefficients:
ar1 ma1
-0.1727 -1.0000
s.e. 0.0308 0.0026
sigma^2 estimated as 106737: log likelihood = -7355.8, aic = 14717.6
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -1.420592 326.547 203.7282 NaN Inf 0.6344327 -0.04113314
Code
summary(west_arima)
Call:
arima(x = west_train_ts, order = c(1, 1, 1))
Coefficients:
ar1 ma1
-0.1698 -1.0000
s.e. 0.0309 0.0026
sigma^2 estimated as 23955: log likelihood = -6593.77, aic = 13193.54
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.6839532 154.6997 101.0808 NaN Inf 0.6424269 -0.02682684
Code
east_values <-forecast(east_arima, h =length(east_test_ts))midwest_values <-forecast(midwest_arima, h =length(midwest_test_ts))south_values <-forecast(south_arima, h =length(south_test_ts))west_values <-forecast(west_arima, h =length(west_test_ts))
### South Regionsouth_initial_arima_acc <- forecast::accuracy(south_values)south_initial_arima_acc %>%kbl(caption ="<center>Initial South Accuracy<center>") %>%kable_styling()
### West Regionwest_initial_arima_acc <- forecast::accuracy(west_values)west_initial_arima_acc %>%kbl(caption ="<center>Initial West Accuracy<center>") %>%kable_styling()
# Grid Search for p and q ------------------------------------------------------p_values <-0:5# AR orderd_values <-0:0# I orderq_values <-0:5# MA order## East - Perform grid search for ARIMA parameters ----------------------------best_east_model <-NULLbest_east_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(east_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_east_aic) { best_east_model <- arima_model best_east_aic <-AIC(arima_model) best_east_order <-c(p, d, q) } } } }}
Code
## Midwest - Perform grid search for ARIMA parameters -------------------------p_values <-0:5# AR orderd_values <-0:1# I orderq_values <-0:5# MA orderbest_midwest_model <-NULLbest_midwest_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(midwest_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_midwest_aic) { best_midwest_model <- arima_model best_midwest_aic <-AIC(arima_model) best_midwest_order <-c(p, d, q) } } } }}
Code
## South - Perform grid search for ARIMA parameters ---------------------------p_values <-0:5# AR orderd_values <-0:1# I orderq_values <-0:5# MA orderbest_south_model <-NULLbest_south_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(south_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_south_aic) { best_south_model <- arima_model best_south_aic <-AIC(arima_model) best_south_order <-c(p, d, q) } } } }}
Code
## West - Perform grid search for ARIMA parameters ----------------------------p_values <-0:5# AR orderd_values <-0:1# I orderq_values <-0:5# MA orderbest_west_model <-NULLbest_west_aic <-Inffor (p in p_values) {for (d in d_values) {for (q in q_values) {# Fit ARIMA model arima_model <-tryCatch( { fit <-Arima(west_ts, order =c(p, d, q)) fit },error =function(e) {NULL } )# Check if ARIMA model was successfully fittedif (!is.null(arima_model)) {# Check AIC criterionif (AIC(arima_model) < best_west_aic) { best_west_model <- arima_model best_west_aic <-AIC(arima_model) best_west_order <-c(p, d, q) } } } }}
### South Regionsouth_tuned_arima_acc <- forecast::accuracy(south_values)south_tuned_arima_acc %>%kbl(caption ="<center>Tuned South Accuracy<center>") %>%kable_styling()
### West Regionwest_tuned_arima_acc <- forecast::accuracy(west_values)west_tuned_arima_acc %>%kbl(caption ="<center>Tuned West Accuracy<center>") %>%kable_styling()
# Building SARIMA model ------------------------------------------------------east_sarima <-arima(east_train_ts, seasonal =list(order =c(1,1,1)))midwest_sarima <-arima(midwest_train_ts, seasonal =list(order =c(1,1,1)))south_sarima <-arima(south_train_ts, seasonal =list(order =c(1,1,1)))west_sarima <-arima(west_train_ts, seasonal =list(order =c(1,1,1)))# Summary of the SARIMA modelsummary(east_sarima)
Call:
arima(x = east_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
0.4706 -0.8016
s.e. 0.0382 0.0211
sigma^2 estimated as 11156: log likelihood = -6206.64, aic = 12419.27
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1606608 105.5715 61.2767 -Inf Inf 0.9398295 0.1140208
Code
summary(midwest_sarima)
Call:
arima(x = midwest_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
-0.3477 -1.0000
s.e. 0.0293 0.0026
sigma^2 estimated as 50173: log likelihood = -6970.99, aic = 13947.98
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.8748943 223.8834 138.9352 NaN Inf 0.6215031 -0.06581815
Code
summary(south_sarima)
Call:
arima(x = south_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
-0.1727 -1.0000
s.e. 0.0308 0.0026
sigma^2 estimated as 106737: log likelihood = -7355.8, aic = 14717.6
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -1.420592 326.547 203.7282 NaN Inf 0.6344327 -0.04113314
Code
summary(west_sarima)
Call:
arima(x = west_train_ts, seasonal = list(order = c(1, 1, 1)))
Coefficients:
sar1 sma1
-0.1698 -1.0000
s.e. 0.0309 0.0026
sigma^2 estimated as 23955: log likelihood = -6593.77, aic = 13193.54
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.6839532 154.6997 101.0808 NaN Inf 0.6424269 -0.02682684
Code
## Forecasting with the SARIMA model -------------------------------------------east_values <-forecast(east_sarima, h =length(east_test_ts))midwest_values <-forecast(midwest_sarima, h =length(midwest_test_ts))south_values <-forecast(south_sarima, h =length(south_test_ts))west_values <-forecast(west_sarima, h =length(west_test_ts))
7.2 Accuracy
Code
## Plotting & Accuracy## East Regioneast_initial_sarima_acc <- forecast::accuracy(east_values)east_initial_sarima_acc %>%kbl(caption ="<center>Initial SARIMA East Accuracy<center>") %>%kable_styling()
## Define the parameter grids --------------------------------------------------p_grid <-0:5# AR parameterd_grid <-0:0# Differencing parameterq_grid <-0:5# MA parameterP_grid <-0:1# Seasonal AR parameterD_grid <-0:1# Seasonal differencing parameterQ_grid <-0:1# Seasonal MA parameters <-12# Seasonal period## East Grid Search ------------------------------------------------------------# Create East Regioneast_train <- east_train %>%select(daily_deaths)# Create an empty data frame to store the resultseast_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(east_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame east_results <-rbind(east_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}d_grid <-0:1# Differencing parameter## Midwest Grid Search ------------------------------------------------------------# Create Midwest Regionmidwest_train <- midwest_train %>%select(daily_deaths)# Create an empty data frame to store the resultsmidwest_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(midwest_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame midwest_results <-rbind(midwest_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}## South Grid Search ------------------------------------------------------------# Create South Regionsouth_train <- south_train %>%select(daily_deaths)# Create an empty data frame to store the resultssouth_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(south_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame south_results <-rbind(south_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}## West Grid Search ------------------------------------------------------------# Create West Regionwest_train <- west_train %>%select(daily_deaths)# Create an empty data frame to store the resultswest_results <-data.frame(order =character(),seasonal =character(),AIC =numeric(),stringsAsFactors =FALSE)# Perform grid searchfor (p in p_grid) {for (d in d_grid) {for (q in q_grid) {for (P in P_grid) {for (D in D_grid) {for (Q in Q_grid) { order <-c(p, d, q) seasonal <-c(P, D, Q, s) model <-tryCatch( { fit <-arima(west_train, order = order, seasonal = seasonal) AIC_val <-AIC(fit)# Append the results to the data frame west_results <-rbind(west_results, data.frame(order =paste(order, collapse =", "),seasonal =paste(seasonal, collapse =", "),AIC = AIC_val)) },error =function(e) NULL ) } } } } }}
# Rebuilding SARIMA Models -----------------------------------------------------east_sarima <-arima(east_train_ts, order =c(5, 0, 4), seasonal =list(order =c(0, 1, 1), period =12))midwest_sarima <-arima(midwest_train_ts, order =c(5, 1, 5),seasonal =list(order =c(0, 0, 1), period =12))south_sarima <-arima(south_train_ts, order =c(5, 1, 5),seasonal =list(order =c(0, 0, 0), period =12))west_sarima <-arima(west_train_ts, order =c(5, 0, 5),seasonal =list(order =c(1, 1, 0), period =12))## Forecasting with the NEW SARIMA model ---------------------------------------east_values <-forecast(east_sarima, h =length(east_test_ts))midwest_values <-forecast(midwest_sarima, h =length(midwest_test_ts))south_values <-forecast(south_sarima, h =length(south_test_ts))west_values <-forecast(west_sarima, h =length(west_test_ts))
Code
## Plotting & Accuracy## East Regioneast_tuned_sarima_acc <- forecast::accuracy(east_values)east_tuned_sarima_acc %>%kbl(caption ="<center>Tuned SARIMA East Accuracy<center>") %>%kable_styling()
# Building AutoARIMA Model -----------------------------------------------------east_arima <-auto.arima(east_train_ts, seasonal =TRUE)midwest_arima <-auto.arima(midwest_train_ts, seasonal =TRUE)south_arima <-auto.arima(south_train_ts, seasonal =TRUE)west_arima <-auto.arima(west_train_ts, seasonal =TRUE)# Summary of the ARIMA modelsummary(east_arima)
Series: east_train_ts
ARIMA(5,1,2)
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2
0.5142 -0.7304 -0.1394 -0.3479 -0.2207 -1.0295 0.7353
s.e. 0.0480 0.0344 0.0426 0.0329 0.0418 0.0410 0.0265
sigma^2 = 6861: log likelihood = -5956.39
AIC=11928.79 AICc=11928.93 BIC=11968.22
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1506172 82.50546 47.1043 NaN Inf 0.7224607 0.01707355
Code
summary(midwest_arima)
Series: midwest_train_ts
ARIMA(5,0,3) with zero mean
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
0.1952 -0.6269 -0.3389 -0.2999 -0.4563 -1.2016 0.9332 -0.3022
s.e. 0.0607 0.0312 0.0356 0.0276 0.0442 0.0702 0.0718 0.0382
sigma^2 = 21712: log likelihood = -6544.65
AIC=13107.31 AICc=13107.49 BIC=13151.67
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.3936445 146.7706 85.18588 NaN Inf 0.3810646 -0.0336952
Code
summary(south_arima)
Series: south_train_ts
ARIMA(5,0,4) with zero mean
Coefficients:
ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
0.1881 -0.9049 -0.0800 -0.4542 -0.5567 -0.9401 1.0118 -0.7340
s.e. 0.0349 0.0292 0.0464 0.0251 0.0303 0.0357 0.0434 0.0525
ma4
0.4897
s.e. 0.0402
sigma^2 = 40663: log likelihood = -6864.82
AIC=13749.64 AICc=13749.86 BIC=13798.93
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1422268 200.7611 119.8731 NaN Inf 0.3732984 -0.06523796
Code
summary(west_arima)
Series: west_train_ts
ARIMA(4,0,3) with zero mean
Coefficients:
ar1 ar2 ar3 ar4 ma1 ma2 ma3
-0.0876 0.1258 -0.5602 -0.4779 -0.5985 -0.3463 0.5343
s.e. 0.0505 0.0450 0.0330 0.0309 0.0509 0.0784 0.0421
sigma^2 = 11991: log likelihood = -6241.45
AIC=12498.9 AICc=12499.04 BIC=12538.33
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.1545496 109.1261 72.50105 NaN Inf 0.4607859 -0.08454213
Code
## Forecasting Test Data -------------------------------------------------------east_values <-forecast(east_arima, h =length(east_test_ts))midwest_values <-forecast(midwest_arima, h =length(midwest_test_ts))south_values <-forecast(south_arima, h =length(south_test_ts))west_values <-forecast(west_arima, h =length(west_test_ts))
8.2 Accuracy
Code
## Plotting & Accuracy## East Regioneast_auto_acc <- forecast::accuracy(east_values)east_auto_acc %>%kbl(caption ="<center>Auto ARIMA East Accuracy<center>") %>%kable_styling()